En este ejercicio vamos a estudiar la tasa de enfermedades transmitidas por los alimentos causadas por el consumo de hortalizas y frutas en la Unión Europea y Estados Unidos con los datos del trabajo de Callejón et al.[1] (2015).

Las enfermedades del análisis fueron Norovirus, Salmonella spp, Escherichia coli, Campylobacter spp, Shigella spp, Clostridium spp, Staphylococcus spp, Yersinia spp, Bacillus spp, otros virus y otros microorganismos. Las verduras eran Ensalada (todos los productos relacionados con la ensalada), Hojas (todos los productos relacionados con las hojas), Tomate y Otras verduras. Las frutas fueron Germinados (todos los productos relacionados con los germinados), Bayas, Melón, Zumos y Otras frutas. Los datos se encuentran en el archivo vegetables.zip. Aunque se pueden cargar desde el archivo vegetables.mat de Matlab, se recomienda leerlos desde los otros archivos. El archivo vegetables.dat contiene tres columnas de datos. La primera columna es la variable illness (enfermedad), la segunda es la variable veg.fruit (verduras/frutas). Y, por último, la tercera columna es la variable region (región). Los nombres de todos los elementos que necesitamos se hallan en el archivo vegetables_labels.txt de forma ordenada según la codificación.

  1. Obtener las tablas de contingencia entre enfermedades y verduras/frutas separadas por región. brevemente.
#cargamos los datos

vegetables <- read.table("vegetables.dat")
noms <- read.table("vegetables_labels.txt")

#eliminamos los datos en blanco
rm.blanks <- function(x) gsub(" ","",x)
noms <- apply(noms,1,rm.blanks)

colnames(vegetables) <- c("illness","veg.fruit","region")
vegetables$illness <- factor(vegetables$illnes, labels=noms[1:11])
vegetables$veg.fruit <- factor(vegetables$veg.fruit, labels=noms[12:20])
vegetables$region <- factor(vegetables$region, labels=noms[21:22])

#hacemos tablas de contingencia para cada region.

attach(vegetables)
tabla.EU <- table(illness[region=="EU"],veg.fruit[region=="EU"])
tabla.EU
##                       
##                        Salad Leafy Tomato Other-vegetables Sprouts Berries
##   Norovirus               15    26      1                9       0      55
##   Salmonella               8    12      1                3      11       0
##   E-coli                   3     0      0                1       3       0
##   Campylobacter            2     1      0                0       0       0
##   Shigella                 1     0      0                3       0       0
##   Clostridium              0     1      0                6       0       0
##   Staphylococcus           0     0      0                3       1       0
##   Yersinia                 0     0      0                3       0       0
##   Bacillus                 2     0      0                3       0       1
##   Other-viruses            2     0      0                0       2       0
##   Other-microorganisms     0     0      0                9       0       0
##                       
##                        Melon Juices Other-fruits
##   Norovirus                0      0            2
##   Salmonella               1      1            3
##   E-coli                   0      0            0
##   Campylobacter            0      0            0
##   Shigella                 0      0            1
##   Clostridium              0      0            0
##   Staphylococcus           0      0            0
##   Yersinia                 0      0            0
##   Bacillus                 0      0            0
##   Other-viruses            0      0            1
##   Other-microorganisms     0      0            0

Tabla de contingencia para USA

tabla.USA <- table(illness[region=="USA"],veg.fruit[region=="USA"])
tabla.USA
##                       
##                        Salad Leafy Tomato Other-vegetables Sprouts Berries
##   Norovirus               97    62      5                9       0       5
##   Salmonella               8     8     17                3      14       2
##   E-coli                  10    22      0                0       4       2
##   Campylobacter            4     2      1                0       0       0
##   Shigella                 1     2      0                0       0       0
##   Clostridium              0     0      0                0       0       0
##   Staphylococcus           2     0      0                0       0       0
##   Yersinia                 0     0      0                0       0       0
##   Bacillus                 1     0      0                0       0       0
##   Other-viruses            0     1      1                1       0       0
##   Other-microorganisms     0     0      0                0       2       0
##                       
##                        Melon Juices Other-fruits
##   Norovirus                9      3           33
##   Salmonella              14      0            5
##   E-coli                   0      6            2
##   Campylobacter            1      0            1
##   Shigella                 0      0            0
##   Clostridium              0      0            0
##   Staphylococcus           0      0            0
##   Yersinia                 0      0            0
##   Bacillus                 0      0            1
##   Other-viruses            0      0            2
##   Other-microorganisms     1      0            0
  1. Realizar un análisis de correspondencias1 de las tablas, representar los resultados y comentarlos brevemente.
library(ca)
ca.EU <- ca(tabla.EU)
round(ca.EU$sv^2, 3) # valores propios
## [1] 0.501 0.425 0.069 0.042 0.035 0.003 0.000 0.000
plot(ca.EU)

Tambien realizamos un plot asimetrico

plot(ca.EU, map = "rowprincipal")

Del mismo modo se puede hacer el otro gráfico asimétrico:

plot(ca.EU, map = "colprincipal")

Realizamos los mismo graficos para USA, pero antes eliminamos datos faltantes

ca.USA <- ca(tabla.USA[-c(6,8),])
round(ca.USA$sv^2, 3) # valores propios
## [1] 0.384 0.139 0.041 0.022 0.008 0.003 0.001 0.000
plot(ca.USA)

Tambien hacemos graficos asimetricos

plot(ca.USA, map = "rowprincipal")

plot(ca.USA, map = "colprincipal")

Para USA los dos primeros componentes explican el 87.2% de los datos, mientras que para la UE explican el 86.1%.

El Norovirus se muestra como responsable de la mayoría de los brotes relacionados con productos agrícolas, seguido de la Salmonela. El Norovirus está relacionado principalmente con el consumo de ensalada en los Estados Unidos y de bayas en la Unión Europea. La Salmonella fue la principal causa de brotes de productos agrícolas en varios estados de Estados Unidos y fue el patógeno implicado en la mayoría de los brotes asociados al consumo de coles, tomates y melones. Como se refleja en los gráficos, el patrón de los brotes de productos frescos difería en Estados Unidos y en la Unión Europea según el tipo de microorganismo y el vehículo alimentario implicado.

Corradino[2] utilizó el MDS para estudiar la estructura de proximidades en una colonia de monos japoneses. Las observaciones se realizaron en un grupo social de 14 monos japoneses durante un período de un año. Los catorce monos se nombran y describen en la Tabla 1. Se observaron las relaciones de proximidad cada 60 segundos. Si dos monos se encontraban a menos de 1.5 metros de distancia y se toleraban, se decía que estaban “cerca”. Se calcularon las disimilaridades para cada par de monos en función del tiempo que la pareja estuvo cerca la una de la otra. Las disimilaridades se estudiaron por separado en la temporada de cría (monk_85.dis) y fuera de ésta (monk_84.dis).

  1. Utilizar las disimilaridades del archivo monk_84.dis2 para definir un objeto del tipo dist en R o una matriz simétrica con esas disimilaridades. Observaremos el orden de la tabla 1 que se puede obtener del archivo monkeys.dat. El vector con las 91 disimilaridades debe rellenar la matriz triangular inferior del objeto por columnas. También habrá que poner nombres a las filas y columnas de la matriz.

Lo mismo para el archivo monk_85.dis con las disimilaridades en la temporada de cría.

El primer paso es crear la matriz de distancia para los datos

monkeys <- read.csv("monkeys.dat", sep=";")
monkeys$age <- factor(monkeys$age)
monkeys$sex <- factor(monkeys$sex)
df <- read.table("monk_84.dis", header = FALSE, sep="", skip=2)
disim <- as.numeric(df$V1[1:91])
mat <- matrix(0, ncol=14, nrow=14)
mat[lower.tri(mat)] <- disim
mat <- mat + t(mat)
colnames(mat) <- rownames(mat) <- monkeys$alias
d.84 <- as.dist(mat)
d.84
##      ALFA FRAN FELL PANC ISA GILD BETI OLGA ORSE ROSS DIVO CIST ELET
## FRAN  126                                                           
## FELL   55  211                                                      
## PANC  367   60   26                                                 
## ISA    55   39   21   12                                            
## GILD   63   39   17   23   1                                        
## BETI   37   74   73   32  11    3                                   
## OLGA   27   33   24    2  28    3    3                              
## ORSE   96   90   68    8  14   26   24    6                         
## ROSS   62   30    8   51  20    1   35    8   42                    
## DIVO   27   18   51  100  10   17    8    5   28    0               
## CIST    4   10  100   38   1    3    3   17   39    3   16          
## ELET    3   19   38   16   0   35    7   16    4    4   62    6     
## EVA    22    5   16   46   8    9    1    6   46   16    9   26  436

Realizamos la matriz de distancias para la otra temporada.

df <- read.table("monk_85.dis", header = FALSE, sep="", skip=2)
disim <- as.numeric(df$V1[1:91])
mat <- matrix(0, ncol=14, nrow=14)
mat[lower.tri(mat)] <- disim
mat <- mat + t(mat)
colnames(mat) <- rownames(mat) <- monkeys$alias
d.85 <- as.dist(mat)
d.85
##      ALFA FRAN FELL PANC ISA GILD BETI OLGA ORSE ROSS DIVO CIST ELET
## FRAN   80                                                           
## FELL   44  161                                                      
## PANC  156  120   27                                                 
## ISA    49   18   10   21                                            
## GILD  109  117   22    4   2                                        
## BETI   77  113   23   40  78   11                                   
## OLGA   83  105    9    9   7   74   18                              
## ORSE   53  153  105   41  25   47   69    4                         
## ROSS   33   44    7   11   4    2   14   16    9                    
## DIVO   20   29   28   75  20   44   60   14   39    0               
## CIST   22   22   20    9   4   56  107   28   62   15   11          
## ELET  286    7   27    7   2   41    7   14    3   15  109   34     
## EVA    59   17  256   18  12   38    0  146   96   20   12   22  204
  1. Comprobar que las disimilaridades en las dos temporadas forman conjuntos no euclídeos. ¿Qué significa esto?
library(ade4)
is.euclid(d.84,print=F)
## Warning in is.euclid(d.84, print = F): Zero distance(s)
## [1] FALSE
#como nos faltan datos chequeamos donde estan

which(d.84==0)
## [1] 54 82

Reemplazamos estos datos por un valor pequeño cercano al "0"

d.84[which(d.84==0)] <- 0.0001
which(d.85==0)
## [1] 70 82
d.85[which(d.85==0)] <- 0.0001

Comprobamos si la distancia es euclidea

is.euclid(d.84,print=F)
## [1] FALSE

En este caso la distancia no es euclidea.

is.euclid(d.85,print=F)
## [1] FALSE

En este caso tampoco es euclidea.

Cuando la disimilaridad es euclídea o una medida directa relativa a ella (tal como el coseno), la mejor elección para el MDS es el MDS métrico. Sin embargo, cuando la disimilaridad no es euclídea o incluso no métrica, debemos decidir entre la solución MDS métrica y la no métrica. El MDS no métrico frecuentemente demuestra ser más realista y mucho mejor en la práctica.

  1. Realizar el MDS más apropiado (clásico, métrico o no métrico) con las disimilaridades de las dos temporadas por separado. ¿Cual es el stress en cada caso?
set.seed(123)
library(MASS)
iso.mds.84 <- isoMDS(d.84)
## initial  value 62.369024 
## final  value 62.351900 
## converged
iso.mds.84$stress
## [1] 62.3519
iso.mds.85 <- isoMDS(d.85)
## initial  value 47.386485 
## final  value 47.345903 
## converged
iso.mds.85$stress
## [1] 47.3459
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
nmds.84 <- metaMDS(d.84, engine="monoMDS", autotransform=FALSE, trace = FALSE)
nmds.84$stress*100
## [1] 20.91301
nmds.85 <- metaMDS(d.85, engine="monoMDS", autotransform=FALSE, trace = FALSE)
nmds.85$stress*100
## [1] 22.84126

Aunque con valores elevados de stress, es evidente que la mejor solución se obtiene con la función metaMDS().

  1. Dibujar la representación de puntos con las mejores soluciones del apartado anterior, una para cada temporada.
plot(nmds.84$points[,1], nmds.84$points[,2], type = "n",
xlab="NMDS1", ylab="NMDS2")
title("Temporada sin cría")
text(nmds.84$points[,1], nmds.84$points[,2],
     col=ifelse(monkeys$sex=="male","blue","red"),
     labels = rownames(nmds.84$points), cex=0.5)

plot(nmds.84$points[,1], nmds.84$points[,2], type = "n",
xlab="NMDS1", ylab="NMDS2")
title("Temporada de cría")
text(nmds.85$points[,1], nmds.85$points[,2],
     col=ifelse(monkeys$sex=="male","blue","red"),
     labels = rownames(nmds.85$points), cex=0.5)

Observamos una mayor proximidad en la representación de la temporada de cría. Los machos aparecen agrupados en las dos configuraciones, aunque con alguna excepción. Las hembras también forman un grupo, es decir, se relacionan más entre ellas en las dos temporadas.

En el trabajo de Corradino se deducía que los machos estaban más cerca en la temporada de cría que fuera de ella. En nuestros gráficos este efecto no aparece.

  1. Comparar las dos configuraciones de las temporadas con la función procrustes().

¿Cual es el mono con menor residuo entre las dos configuraciones?

proc <- procrustes(nmds.84,nmds.85)
plot(proc)

El menor residuo de la representación es para Gilda.

sort(residuals(proc))[1:3]
##     GILD      EVA     DIVO 
## 31.67155 41.72386 61.20314
  1. Como el stress es alto, podemos optar por una representación en tres dimensiones.

En primer lugar procedemos a calcular el MDS no métrico con 3 coordenadas.

nmds.85.3 <- metaMDS(d.85, engine="monoMDS", k = 3, autotransform=FALSE, trace = FALSE)
nmds.85.3$stress*100
## [1] 14.81673

El stress es más bajo que con dos dimensiones.

library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
datos <- data.frame(nmds.85.3$points)
fig3D <- plot_ly(data=datos, x = ~MDS1, y = ~MDS2, z = ~MDS3, 
                 color = ~monkeys$sex, colors = c("red","blue"))
fig3D
## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Con la rotación adecuada se puede ver la separación total entre machos y hembras.

Normalmente, el análisis de componentes principales se lleva a cabo calculando los valores propios y los vectores propios de la matriz de correlaciones. Con n casos y p variables, si escribimos X para la matriz n×p estandarizada para que las columnas tengan media cero y desviación estándar unitaria, encontramos los valores propios y los vectores propios de la matriz p × p X0X (que es n o n − 1 veces la matriz de correlaciones dependiendo de qué denominador se haya utilizado al calcular las desviaciones estándar). Con los datos de los gorriones del ejercicio 3 de los ejercicios propuestos para PCA, comprobar las siguientes afirmaciones:

  1. El primer vector propio da las cargas (loadings) de cada variable en la primera componente principal, el segundo vector propio da las cargas en la segunda componente y así sucesivamente. Comprobar este resultado con la función princomp().
load("gorriones.RData")
X <- gorriones[,1:5]
X <- scale(X)
eig <- eigen(crossprod(X), symmetric = TRUE)
pca <- princomp(X)
pca$loadings
## 
## Loadings:
##    Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## x1  0.452         0.690  0.420  0.374
## x2  0.462 -0.300  0.341 -0.548 -0.530
## x3  0.451 -0.325 -0.454  0.606 -0.343
## x4  0.471 -0.185 -0.411 -0.388  0.652
## x5  0.398  0.876 -0.178        -0.192
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## SS loadings       1.0    1.0    1.0    1.0    1.0
## Proportion Var    0.2    0.2    0.2    0.2    0.2
## Cumulative Var    0.2    0.4    0.6    0.8    1.0
  eig$vectors
##           [,1]        [,2]       [,3]        [,4]       [,5]
## [1,] 0.4517989  0.05072137  0.6904702  0.42041399 -0.3739091
## [2,] 0.4616809 -0.29956355  0.3405484 -0.54786307  0.5300805
## [3,] 0.4505416 -0.32457242 -0.4544927  0.60629605  0.3427923
## [4,] 0.4707389 -0.18468403 -0.4109350 -0.38827811 -0.6516665
## [5,] 0.3976754  0.87648935 -0.1784558 -0.06887199  0.1924341
  1. Escribiendo las cargas de las primeras q componentes como columnas de la matriz p×q B, la matriz n × q P de las puntuaciones (scores) de las componentes principales de los individuos se obtiene aplicando las cargas a la matriz de datos original, es decir, P = XB.
scores <- X %*% eig$vectors
head(scores)
##             [,1]        [,2]       [,3]         [,4]       [,5]
## [1,]  0.06428901 -0.60083713 -0.1712334 -0.515825561  0.5487904
## [2,] -2.18031283 -0.44230082  0.4000696 -0.645459959  0.2310766
## [3,] -1.14556567  0.01925412 -0.6761269 -0.716298164  0.2088714
## [4,] -2.31106565  0.17199267 -0.3059621  0.149289289  0.4781034
## [5,] -0.29504203 -0.66520783 -0.4742138 -0.545862110  0.2444780
## [6,]  1.91626198 -0.59525444  0.6209330  0.006608669 -0.2855166
head(pca$scores)
##           Comp.1      Comp.2     Comp.3       Comp.4     Comp.5
## [1,]  0.06428901 -0.60083713 -0.1712334 -0.515825561 -0.5487904
## [2,] -2.18031283 -0.44230082  0.4000696 -0.645459959 -0.2310766
## [3,] -1.14556567  0.01925412 -0.6761269 -0.716298164 -0.2088714
## [4,] -2.31106565  0.17199267 -0.3059621  0.149289289 -0.4781034
## [5,] -0.29504203 -0.66520783 -0.4742138 -0.545862110 -0.2444780
## [6,]  1.91626198 -0.59525444  0.6209330  0.006608669  0.2855166
  1. La matriz de suma de cuadrados y productos, P'P = D, es diagonal con elementos iguales a los primeros q valores propios de de la matriz X0X, de modo que las varianzas de las componentes principales pueden obtenerse dividiendo los valores propios por n o n − 1.
D <- crossprod(scores)
eig$values
## [1] 173.566961  25.512196  18.548378  14.475145   7.897321
diag(D)
## [1] 173.566961  25.512196  18.548378  14.475145   7.897321
diag(D)/49
## [1] 3.5421829 0.5206571 0.3785383 0.2954111 0.1611698
pca$sdev^2
##    Comp.1    Comp.2    Comp.3    Comp.4    Comp.5 
## 3.5421829 0.5206571 0.3785383 0.2954111 0.1611698
  1. Sustituimos los valores desconocidos (missings) de la matriz X (la de los controles) por la media de cada columna. El siguiente es estandarizar los datos para que tengan media cero y desviación estándar uno y finalmente calcular la matriz XX'.
load("SNP.RData")
controls <- rownames(subject.support)[subject.support$cc==0]
pop <- subject.support[controls,"stratum"]
pop.all <- subject.support[,"stratum"]
use <- seq(1, ncol(snps.10), 10)
## Loading required package: snpStats
## Loading required package: survival
## Loading required package: Matrix
ctl10 <- snps.10[controls, use]
all10 <- snps.10[,use]
X <- as(ctl10, "numeric")
X.all <- as(all10,"numeric")
dim(X)
## [1]  500 2851

La primera dificultad es substituir los valores missing por la media.

for(i in 1:ncol(X)){
  X[is.na(X[,i]), i] <- mean(X[,i], na.rm = TRUE)
}

Estandarizamos los datos

X <- scale(X)

Calculamos la matriz XX'

XXt <- tcrossprod(X)
head(XXt[1:10,1:10])
##           jpt.869   jpt.862   jpt.948   ceu.564   ceu.904   ceu.665   jpt.663
## jpt.869 2601.7764  138.9807  192.9046 -283.5339 -414.6340 -235.4016  247.0385
## jpt.862  138.9807 2775.1951  164.8674 -218.8962 -306.4598 -280.7771  202.4891
## jpt.948  192.9046  164.8674 2531.5733 -361.6866 -311.3615 -278.4234  309.6210
## ceu.564 -283.5339 -218.8962 -361.6866 2908.3334  439.0646  329.4606 -162.9436
## ceu.904 -414.6340 -306.4598 -311.3615  439.0646 2935.7316  313.8808 -296.6975
## ceu.665 -235.4016 -280.7771 -278.4234  329.4606  313.8808 2836.9213 -326.4090
##           ceu.977   jpt.637    ceu.897
## jpt.869 -311.0834  329.9604  -95.97765
## jpt.862 -345.0864  304.8939 -306.84784
## jpt.948 -284.2715  316.8407 -261.73112
## ceu.564  260.8681 -310.9871  160.09227
## ceu.904  229.0223 -403.5789  308.65566
## ceu.665  355.8414 -296.9524  329.36431
  1. Calcular los valores y vectores propios de la matriz XX' y, a partir de ellos, la matriz U. Comparar gráficamente las puntuaciones de la primera componente para los controles según la población.

Hacer el mismo gráfico con la segunda y tercera componentes. ¿Qué componente separa mejor las poblaciones?

eig <- eigen(XXt, symmetric = TRUE)
eig$values[1:5]
## [1] 146789.918   8589.204   8380.229   8236.379   8111.603
U <- eig$vectors
oldpar <- par(mfrow=c(1,3))
boxplot(U[,1] ~ pop, ylab="PC1")
boxplot(U[,2] ~ pop, ylab="PC2")
boxplot(U[,3] ~ pop, ylab="PC3")

par(oldpar)
  1. Calcular las cargas B según se ha explicado antes.
eig$values[eig$values<0] <- 0
v <- 1/sqrt(eig$values)
m <- crossprod(X, U)
# B <- m %*% diag(v) # poco eficiente
# B <- t(t(m) * v) # más eficiente
B <- m * rep(v, rep(nrow(m), ncol(m))) # mejor
B[1:10,1:10]
##                     [,1]          [,2]         [,3]         [,4]         [,5]
## rs7909677   0.0003252226 -9.331756e-03  0.015445805  0.003947906 -0.011546724
## rs4880781  -0.0055876444  7.507270e-03 -0.000543355 -0.008346595 -0.003828963
## rs2448365  -0.0076646829  1.733169e-02 -0.008138940  0.007698970  0.007641585
## rs7919436   0.0209014311  1.102595e-02 -0.025856734  0.013353554  0.024974745
## rs12357593  0.0301087215  1.708045e-03 -0.022693905  0.014494197  0.022537829
## rs10904173 -0.0032615624  1.569638e-02  0.017220464  0.002238865 -0.014006605
## rs11252693  0.0107882727 -4.632472e-05 -0.012653201  0.011259890  0.006176045
## rs11253096  0.0259380458 -1.762775e-02  0.006504495  0.010741717 -0.010751181
## rs17159711  0.0150469759  1.284841e-02  0.010802802  0.007726009  0.030637888
## rs4881518   0.0168951870 -1.775938e-02 -0.008888053  0.018592907  0.004927426
##                    [,6]         [,7]          [,8]         [,9]        [,10]
## rs7909677   0.009234650  0.009592387 -0.0006004881  0.009523450 -0.014859529
## rs4880781  -0.025199111  0.001262191  0.0136031063 -0.034072213 -0.011145068
## rs2448365  -0.009577234  0.025332822  0.0176713475  0.005336778 -0.009796594
## rs7919436   0.004708050 -0.017783670  0.0129060916  0.035684897 -0.009938182
## rs12357593  0.017645674 -0.035232563  0.0028693908  0.036582396 -0.002629847
## rs10904173  0.008486269  0.013840947  0.0097519493 -0.023745230 -0.017696460
## rs11252693  0.019020146 -0.019180868 -0.0004079738  0.033007216  0.014243545
## rs11253096  0.004788495 -0.014961988  0.0022505070 -0.024822305 -0.008882354
## rs17159711 -0.011462165  0.022790028 -0.0054857879  0.008332284 -0.009762654
## rs4881518   0.022255034 -0.006821502  0.0213503830  0.025230784  0.001849588
  1. Calcular las puntuaciones de las dos primeras componentes de todas las observaciones con las cargas de los controles.

Se controlan los valores faltantes y se estandarizan los datos

for(i in 1:ncol(X.all)){
  X.all[is.na(X.all[,i]), i] <- mean(X.all[,i], na.rm = TRUE)
}
X.all <- scale(X.all)

Calculamos las puntuaciones y hacemos el grafico

scores.all.12 <- X.all %*% B[,1:2]
plot(scores.all.12, xlab="PC1", ylab="PC2",
     col=ifelse(pop.all=="CEU","red","blue"))

Las dos poblaciones quedan también separadas por la primera componente.